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Abstract 

In this paper numerical methods of computing distances between two Radon mea- 
sures on K are discussed. Efficient algorithms for Wasserstein-type metrics are pro- 
vided. In particular, we propose a novel algorithm to compute the flat metric (bounded 
Lipschitz distance) with a computational cost 0(n log n). The flat distance has recently 
proven to be adequate for the Escalator Boxcar Train (EBT) method for solving trans- 
port equations with growth terms. Therefore, finding efficient numerical algorithms to 
compute the flat distance between two measures is important for finding the residual 

£T} error and validating empirical convergence of different methods. 

t— I Keywords: metric spaces, flat metric, Wasserstein distance, Radon measures, 

optimal transport, linear programming, minimum-cost flow 
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1 Introduction 



Recent years witnessed large developments in the kinetic theory methods applied to mathe- 
matical physics and more recently also to mathematical biology. Among important branches 
of the kinetic theory are optimal transportation problems and related to them Wasserstein 
metrics, and Monge-Kantorovich metrics [21 [Hj. Partial differential equations in metric 
spaces are being applied to transportation problems |10| 114] . gradient flows [2j [13] and 
structured population models (H [7J El I12| . 
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Output of mathematical modeling can often be described as Radon measures. Compar- 
ing the results of the models requires then a definition of distance in the space of measures. 
The desired properties of such metrics depend on the structure of the considered problem. 
In most of the cases, the topology of total variation is too strong for applications, and 
weaker metrics had to be introduced, see [8] for details. Well known 1-Wasserstein met- 
ric, on the other hand, is only applicable to processes with mass conservation. To cope 
with growth in the process, various modifications have been proposed, including flat met- 
ric and centralized Wasserstein metric. In the present paper, we additionally introduce a 
normalized Wasserstein distance. For comparison of different metrics, their interpretation 
and examples see Table [T] Even though, all those distances can be computed using linear 
programming (LP), its computational complexity becomes often larger than the complex- 
ity of solving the original problem. For example, the equations for which the stability of 
numerical algorithm is proven in (W 1,co )* require an efficient algorithm for the flat metric 
to find the residual error. 

Algorithms proposed in this paper are designed to compute efficiently Wasserstein-type 
distance between two Radon measures in the form of Y17=i a i$xii where 5 is the Dirac delta. 
The algorithms can be also applied for the case of an arbitrary pair of measures, by ap- 
proximating those measures by the sum of Dirac deltas. In the process of finding numerical 
solutions to partial differential equations, the distance between a discrete and absolutely 
continuous measures is needed to evaluate the quality of the initial condition approximation, 
while the distance between two discrete measure is needed to find the residual error. 

The paper is organized as follows. In Section [2] we introduce and compare four Wasserstein- 
type metrics. Additionally, we introduce a two-argument function, based on the Wasserstein 
distance, which is not a metric, however, provides a good tool to estimate the flat metric 
from above. Section [3] is devoted to numerical algorithms proposed to calculate distances 
between two measures in respect to the considered metrics. We justify the algorithms and 
provide respective pseudocodes. The novelty of this paper is the algorithm for flat met- 
ric, which has been recently proven to be adequate for the Escalator Boxcar Train (EBT) 
method for solving transport equations [3]. We propose an efficient algorithm which com- 
putation cost is 0(n 2 ) and optimize it further to 0(n log n). To judge the efficiency of the 
algorithms, we compare the times needed to compute the flat distance between two sums 
of a large number of Dirac deltas randomly distributed over [-1,1]. 

2 Metrics on the spaces of Radon measures 

2.1 1- Wasserstein distance 

The framework of Wasserstein metric in the spaces of probability measures has proven 
to be very useful for the analysis of equations given in a conservative form, such as, for 
example, transport equation |2j, Fokker-Planck [9] or nonlinear diffusion equation [5]. It 
was originally defined using the notion of optimal transportation, but we focus on its dual 



representation. 

Definition 2.1. 1-Wasserstein distance between two measures \x and v is given by 



Wx(fj,,u) =sup 



f(x)d{ii-v){x) 



:/eC(l,l),L f (/)a 



Observe that the integral J R f(x)d(fi — v){x) is invariant with respect to adding a constant 
to function /. Consequently, for any i*GM holds 

Wxfov) = supjy f(x)d(n-v)(x) : / G C(R,R),/(ar*) = 0,Lip(f) < l| . 

Furthermore, the following two properties hold: 

• 1-Wasserstein distance is scale-invariant 

Wi(X ■ fi, X ■ v) = AWi(/i,i/). 

• 1-Wasserstein distance is translation-invariant 

2.2 Normalized Wasserstein metric 

If fJ,(M) 7^ ^(R), then Wi(fi,u) = oo by the dual representation definition. It makes this 
metric useless for processes not preserving the conservation of mass. 

In this section we introduce a simple normalization that leads to a definition of translation- 
invariant metric suitable for non-local models 

Definition 2.2. We define normalized 1-Wasserstein distance between two measures /i and 
v as 

Wiin, u) = min ( H/xll + H , | IHI " INI I + Wi (^, ~X) , (2) 

where \\fi\\ denotes the total variation of the measure fj>. 
Lemma 2.3. The distance defined by (pi) is a metric. 
Proof. Let fi, v and 77 be Radon measures. Then, it holds 
1. W\(n, v) = if and only if jjl = v. 

Indeed, either ||//|| + \\v\\ = or I llull — IMI I + W\ ( w^r, w^r ) = imply that u = v. 



2. W 1 (ji,v) = W 1 (v,n), 

3. Since 

Wi(fj,,v) + Wi(u,rj) = min ( \\/j,\\ + \\v\\ , | ||/x|| - \\i/\\ | + W\ ( Tr - 7, 7— r 

V VIImII IM 

/.. , ,, / Tj V 

+ min M + M>M — M + Wi 71 — it > 71 — r 
to show the triangle inequality, we consider four possibilities 
Wi(/x, 1/) + Wi(v,r]) = \\fi\\ + ||i>|| + ||r/|| + ||i/|| ^ ||/i|| + ||r/|| ^ Wi(//,77), 

WiG^i/J + Wi^) = HlMll-lklll+^ifA^o 

VIlMll IM 

+ 1 Ihll - H I + W! (-^r ^r) > W^n), 

\\\v\\ HI/ 

Wi(p,u) + Wi(v,T)) = \\fi\\ + \\v\\ + I ||r?|| - \\v\\\ + Wi (7-^77,77^77 ) ^ ||m|| + IN 

VIMI HI/ 

> Wi(M,r?), 

W 1 {jl,u) + Wi{u,r]) = \\f]\\ + ||l/|| + | ||/x|| - ||l/||| +Wi ( ipTT, tt^tt ) > \\fi\\ + \\7)\ 

VIImII IMI/ 

> Wi(p,i/). 



□ 

This metric can be easily computed numerically using the algorithm for 1-Wasserstein 
distance. However, it lacks the scaling property, which holds for the Wasserstein distance, 
and which is useful for applications. Nonetheless, the following weaker property holds. 

Proposition 2.4. Let [x k and v k be two sequences of Radon measures and \\/J-k\\ -^ 0, 
\\i/ k \\ ->■ then Wi{n kl v k ) ->■ 

2.2.1 Centralized Wasserstein metric 

In this section we present a different modification of the 1- Wasserstein distance, which is 
scale-invariant. 

Definition 2.5. We define centralized 1-Wasserstein distance between two measures /i and 
v as 



Wi(fi, v) = sup 



f(x)d(fi- v){x) 



:/€C(R,R),Lip(/)<l,/(0)€[-l,l]L 



The metric was introduced in [H] for analysis of the structured population models in 
the spaces of non-negative Radon measures, A^^M" 1 "), on M + with integrable first moment, 
i.e., 

M!(R + ) := {n e M + (R + ) I f \x\ dfj, < oo). 
L Jr+ j 

This metric satisfies the scaling property, but on the contrary to Wasserstein metric, it 

is not invariant with respect to translations. It is therefore only useful for modeling specific 

phenomena, for example structured population dynamics where mass generation can only 

occur in one specific point of space. 

2.3 Flat metric 

Another solution of the problem of comparing two measures of different masses is the flat 
metric which is defined as follows 

Definition 2.6. Flat distance between two measures \i and v is given by 



F(/u, u) = sup 



f{x)d{n-v){x) 



:feB c{%m) (0,l),Lip(f)^l 



The flat metric, known also as a bounded Lipschitz distance [11], corresponds to the 
dual norm of W 1,00 (]R), since the test functions used in the above definition are dense 
in M /1,oc (]R). This metric satisfies the scaling property and it is translation invariant. It 
has proven to be useful in analysis of structured population models and in particularly, 
Lipschitz dependence of solutions on the model parameters and initial data [2 HI- The flat 
metric has been recently used for the proof of convergence and stability of EBT, which is 
a numerical algorithm based on particle method, for the transport equation with growth 
terms [3]. Consequently, an implementation of the EBT method requires an algorithm to 
compute the flat distance between two measures. 

2.4 Upper bound for flat metric 

Computing of flat metric exactly is more costly than computing of 1- Wasserstein metric 



(0(n log n) vs 0(n), see Section 3.2.4). Moreover, often in applications it is sufficient to 
calculate upper bound of the distance, for example when computing residual error. We 
propose the following function, which requires only linear computing time. 



F(ji,v) = \\\i4 -\M\ + { > N 



if ||/x|| < \\u\ 
if \\ji\\ > \\v\ 



Lemma 2.7. Let \x, v be Radon measures on R. Then, the following estimate holds 

F{frv)<F{n,v) 



Proof. Assume that 1 1 ^11 < ||zv||. Then, using the estimate of the flat metric by 1-Wasserstein 
distance for the measures with equal masses provides 



F(ji,v) < F(ji,v 



\fi 



= wAn.v 



+ F(i/-Mi/) 



+ F(u 
+ sup 

+ 



I Ml I 



/(*)(■ 



w 



l)dv(x) 



:/eB CW) (0,l),If(/)O 



IMII 

Ml 



The last equality results from the representation of Radon distance and positivity of measure 
v. □ 

The upper bound function is more useful if it can be estimated from above by c-F(/j,, v) 
for some constant c. In a general case, however, such constant does not exist, since by 
taking \i = Sq, v = 5 X and passing x — > oo we obtain 



F(jjt, V ) 
F(ji,u) 



oo, 
2. 



Nevertheless, the desired estimate can be shown on a bounded set. 

Lemma 2.8. Let /j,, v be non-negative Radon measures on a compact set K, then the 
following estimate holds 

F{n,v) < c K F{n,v) 

Proof. Let ck = max(l, ^diam^K)) and xq be such point that \k — Xq\ < ck for any k £ K. 
Then, assuming ||/x|| < ||i/||, we obtain 



F{n,v) 



w\ 



\u\\ \+Wl[lJL,V 



I Ml 



r d ( P -v)\ +Wl L^\ 



< F^, v) + W 1 [ix,v 



I Mil 



< F^, v) + W x (m, v) + W 1 [v,v 



I Ml 



< 2-F{ l i,v) + W 1 {ii,v) 
= 2 • F(ji, v) + sup 

= 2 • F(fi, v) + c_ft: sup 

< (2 + c^)F(/x, i/) 



f(x)d(n-u)(x) 



:/eC(R,R),Lip(/)<l 



f{x)d{n-v)(x) 



:f£C(R,R),Lip(f)^ 



c K 



□ 



2.5 Comparison of the presented metrics 

In the following Table we compare the introduced metrics and provide examples and inter- 
pretations in terms of optimal transport. 



Metric 


Example: 


Scale- 


Translation- 


Intuition of d(fi, v) 


Compute 




d(25 x ,35 y ) 


invariancc 


; invar iance 




complexity 


Wasser- 


oo 


YES 


YES 


The cost of optimal 


0(n) 


stein 








transportation of distribution 

jutoa state given by v, 
assuming that moving mass 
m by x requires rax energy. 





Wasser- 
stein 

normal- 
ized 


min(2 + 3, 
(3-2) + |x-y|) 


weak 


YES 


Minimum of the cost of 

annihilating [i and generating 

v; and of the cost of 

generating the difference in 

mass between \i and v and 

transporting #tt to yr+, 

assuming that 

generating/annihilating mass 

at any position requires the 

cost equal to the mass. 


0(n) 


Wasser- 
stein 

central- 
ized 


2\x-y\ + \y\ 


YES 


NO 


The cost of generating the 

difference in mass in point 

in space added to the cost of 

transporting 

M+ (INI - IImID^o to v. 


0(n) 


Flat 


1 + 
2min (2, \x — y\) 


YES 


YES 


The cost of optimal 

transporting AND/OR 

generating AND/OR 

annihilating mass to form v 

from ji. 


0(n\ogn) 


Radon 


2 + 3 


YES 


YES 


The cost of generating 

AND/OR annihilating mass 

to form v from fi 


<D(n) 






!able 1: Co 


mparison of c 


ifferent metrics. 





3 Approximations and computational algorithms 

In the remainder of this paper, we present algorithms for computing and approximating of 
the introduced metrics. 



3.1 Computing of 1-Wasserstein distance 

We start with introducing the algorithm for computation of 1-Wasserstein distance. The 
algorithm is well-known and straightforward. Nevertheless, we present it here, since the 
specific approach used in this section will be also applied later in the proof of correctness 
of more involved algorithms for other distances. 



3.1.1 Reduction to the case of discrete measures 

Lemma 3.1. Let [i be an arbitrary Radon measure supported on [0,1]. For every e > 
there exist a discrete measure [i £ such that 



Wi(n,ii e ) <e 



Proof. Define 



/ Mo,i] 



V^ 8i/nV 



1=1 



i — 1 i 



n n 



Then, we estimate 






i=l 



i — 1 i \ 1 r 

,- <-m0,1, 

n n I n 



and the right-hand side tends to with n — > oo. 



□ 



Consequently, 1-Wasserstein distance between an arbitrary pair of finite Radon measures 
can be estimated by the distances between their discrete approximations. From this point 
on, we assume that 



fc=i 



3.1.2 The algorithm 

The formula for 1-Wasserstein distance reads 



W^u) = sup \ Y^a k f{x k ) : / € C(R,R),f(x n ) = 0,Lip(f) < 1 I . 



.fc=i 



Regularity conditions can be represented as linear programming bounds. Hence, computing 
of Wi(n, v) is equivalent to finding maximum of 



y2 a kfk 



fc=i 



with the following restrictions 



I fk ~ fi 



fc-i 



< 



l-^fc 3?fc — 1 1 



Although this problem can be solved by linear programming, a more efficient algorithm can 
be found. 
Define 

W^{x) = SUp I ^ a kfk ■ fm= X, \f k - / fe _i| ^ \x k - X k -l\ \ . 

Then, obviously Wi(p, v) = VF"(0). Denote d k = Xfc+i — x k , and observe that 

Wl(x) = a\x, 

Wf(x) = a>2X + sup Wl(x) = a>2X + a\x + ai- sgn(ai)d\ = 

fie[x—di,x+di] 

= (ai + a 2 )x + \ai\d\. 



It can be shown by induction that 

/ n \ n—l 






1 

J2 a i 
3=1 



Notice that the value a n is not used in this formula. It is, however, involved indirectly, as 
EIli a t = 0. 

3.1.3 Pseudocode 

Input: Non-decreasing table of positions x£ 
[0, l] n , table of masses a £ W l 
1-Wasserstein-Distance (x £ [0, l] n , a G R n ) 

distance ^— 

partialSum ■<— 

for idx*r- 1 to n — 1 do 

partialSum ■<— partialSum + a^ x 
distance <— distance + (xjdx+i — x^x) • \partialSum\ 
return distance 

3.1.4 Complexity of the algorithm 

It is clear from the pseudocode that the computational complexity of the algorithm is B(n), 
while memory complexity is 0(1). 
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3.2 Computing of the centralized Wasserstein distance 

3.2.1 Reduction to the case of discrete measures 

Similarly to the case of 1- Wasserstein distance the following lemma holds. 

Lemma 3.2. Let fj, be an arbitrary Radon measure supported on [0,K]. For every e > 0, 
there exists a discrete measure jjl £ such that 

Wl(fJL,(l E ) < £. 

The proof follows from the fact that ||//|| = ||^|| implies W\(fj,, v) = Wi(/x, v). 

3.2.2 The Algorithm 

We assume 



A* - v = Z ai ^ x i + a ™+i^o + Z a,i6 Xi . 

i=l i=m+2 



Define 



W{(x) = sup I Y^ o-kfk ■ fj = x, \f k - fk-i\ < 



%k x k — 1| i > 



,k=i 



(Fj(x) - sup { ^ a kfk ■ fj = x,\fk- fk-l\ < \xk - x k -i\ 

k=j 



As already proven 



fm-\-\ 



Wf +i (x) 



W 1 (x) 



fc=i 



y^ aj ) x + y^ d k 

n 

E 



E 

i=\ 



i=l / 

n \ n— (m+1) 



\i=m+l 

From LP representation of the metric 



fc=i 



i=rz+l— k 



Wi(fi,u) = sup <^ ^a k f k : -1 < / m+i < 1, \f k - / fe _i| < |ar fe - x fc -i| 



,fc=i 



it can be deduced that 



W 1 { f i,v)= sup [W™ + \x) + W'r (x)-a rn+1 x), 
xe[-i,i] 



11 



so the distance is given by the formula 



m 


k 


n— (m-fl) 


n 




n 


iO,^) = Y dk 


Y Qi 


+ 2^ ^n-fc 


S °' 


+ 


J2 Qi 


fe=l 


i=l 


fc=l 


i=n+l— k 




i=l 



3.2.3 Pseudocode 



Input: Non-decreasing table of positions i£ 



-fc-l 



table of masses a € 



[0,lf x{0} x [0,1]"' 
1-Wasserstein-Central i zed-Distance 

distance <— 

(partial SumFront, partial SumBack) <— (0, 0) 

(idxFront, idxBack) <— (l,n) 

while x idxFront < do 

partialSumFront <— partial SumFront + aidxFront 
distance <- distance + (xidxFront+i - XidxFront) ■ 
\partialSumFront\ 
idxFront <— idxFront + 1 

while XidxBack > do 

partial SumBack <— partial SumBack + aidxEnd 
distance <- distance + (x idx Back ~ XidxBack-i) ■ 
\partialSumBack\ 
idxBack ^— idxBack — 1 

for idx •(— idxFront to idxBack do 

partialSumFront ■<— partialSumFront + a^x 
return distance + {partial SumFront + partial SumBack\ 



3.2.4 Complexity of the algorithm 

Each iteration of each loop takes a constant time. The total number of iterations in all 
three loops is equal to k + \ + (n — k — 1). Computational complexity of this algorithm is 
therefore 0(re), while the memory complexity is 0(1). 



3.3 Computing of flat distance 

The algorithm for flat distance requires storing the shape of functions analogous to W" 



as they get more complicated when m increases. In Section 3.3.2 we provide a recursive 
formula for the sequence of these functions. The pseudocode in Section |3.3.3 implements 
the algorithm using an abstract data structure to store previously defined functions. The 
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computational complexity depends on the choice of this structure. In further sections we 
provide two solutions that require 0(n 2 ) and O(nlogn) operations. 

3.3.1 Reduction to the case of discrete measures 

Similarly to the previous cases the following lemma holds. 

Lemma 3.3. Let \i be an arbitrary Radon measure supported on [0, K\. For every e > 
there exist a discrete measure [i £ such that 

F(ji,fJte) <£ 

The proof follows from the fact that ||/x|| = \\v\\ implies F(fi, v) < Wi([i, u). 

3.3.2 The Algorithm 

We assume 

n 
[I- V = y^ajJa;-. 
i=l 

Computing of F(p, v) is equivalent to finding maximum of 

n 

/ , Qfc/fc 



fc=l 



with the following restrictions 



l/fcl < 1, 

\fk-fk-i\ < \Xk~Xk-l\- 



Define 



F m {x) = sup 



yZ ak f k 



k=i 



■ fm = X, \fk ~ fk-l\ < \Xk ~ Xk-l\ , \fk\ < 1 ) ■ 



Obviously it holds by definition 



F(/jt,v)= sup F n (x) 

x6[-l,l] 
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Observe that 

F 1 (x) = a\x, 

F 2 (x) = (12X + sup F 1 (x) = a,2X + min(|ai|, a\x + |ai|di), 

/lS[x-di,x+di]n[-l,l] 



F m (x) = a m x+ sup F m - 1 (x). 

/ m _i€[x-d m _i,x+d m _i]n[-l,l] 

Computing of F m (x) based on F m ~ l is more complex than in previous cases, as F m ~ 1 is 
not necessarily monotonic. 

Lemma 3.4. Function F m is concave for each m. 

Proof. To prove the lemma we will use induction with respect to m. F 1 (x) is given as a±x, 
so it is indeed concave. Assume F m is concave. Define 

iCi(x) = sup F n (x). 

[x-d,x+d]n[-l,l] 

Choose x,y G [—1,1]. Then, there exist x' G B(x,d) fl [—1,1], y' G B(y,d) D [—1,1] such 
that 

«F™i(x) + (1 - a)F^(y) = «F™(,') + (1 - a)i^(y')- 
Because F" 1 is concave, it holds 

aF m {x') + (1 - a)F m (y') < F m (ax' + (1 - a)y') ^ F™±(ax + (1 - o)y) 

The last inequality follows from ax 1 + (1 — a)y' G B(ax + (1 — a)y, (i). It is now proven that 
i ?m+1 (x) is convex, as it is a sum of a linear function and a convex function Fmax(y)- O 

Lemma 3.5. For each m the function F m is piecewise linear on m intervals and it holds 
for some point x m 



a m x+ < 



F m 1 (x + d m -i) on[-l,x m -d m -i] 

■F \%m) On \X rn UfYi—l, x rn -\- U m — ij 

F m - l ( x -d m -i) on [x m + dm-i , 1] 



Proof. F l is a linear function, so it can be described by its values in ±1. Assume that F m 
can be described by at most m + 1 points and is linear between those points. As F m is 
concave, there exists a point x m G [— 1, 1] such that F m (x) < F m (x m ) for every x. The 
maximum of F m on an interval whose both ends are smaller than x m is taken on its right 
end. Similarly, if both ends of the intervals are larger than x m , the maximum is taken on its 
left end. Finally, if the interval contains x m , the maximum is in x m . These considerations 
prove the formula for F m+1 . It follows that F m+ \ is piecewise linear and it can be described 
by as many points as F n \, d 1 _ d -> plus 1. □ 
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3.3.3 Pseudocode 

In the following code a data structure called 'funcDescription' being a set of pairs will be 
used to describe F tdx . Its has the following interpretation: 

1. F idx (-l) = leftValue 

2. if (v,p) £ funcDescription then F idx (x)' = p for all x larger than v and smaller than 
the next value in the structure. 

The last instruction in the main loop of the pseudocode makes it inefficient to implement 
'funcDescription' as a simple BST tree. 

Input: Non-decreasing table of positions x£ 
[0, l] n , table of masses a £ 1" 
Flat-Distance 

leftValue <— 

funcDescription <— {(—1, 0), (1, — oo)} 

for idx <- 1 to n do 

& > %idx %idx—l 

funcLeft <— 

{(v — d,p) : (v,p) £ funcDescription f\p > 0} 

funcRight «— 

{(v + d,p) : (v,p) £ funcDescription Ap < 0} 

v m <— minju : (v,p) £ funcRight} 

funcDescription <— funcLeftL){(v m — 2d, 0)}U funcRight 

leftValue <— leftValue + 

E(« lP )e/«ncDe a cripMof»,t;<i (min(min{*/ : (t/, _) £ funcDescription Av' >v}, -1) - «)• 

P 

\VminiPmin) > 

max{(f,p) : (f,p) £ funcDescription[i] Av ^ —1} 

\VmaxiPmax) * 

max{(t>,p) : (f,p) £ funcDescription[i] Av £ [0,1]} 

funcDescription <— 

(funcDescription n {(w,p) : -u £ [0, 1]}) U 

{(max(-u m i„, -l),Pmm)} U {(1, -oo)} 

funcDescription ■<— 

{(x,p + a„) : (x,p) £ funcDescription} 

return leftValue + 

E(^p)e/ u „ c Des C ription,p>o ( min K : K> -) G funcDescription A v' > v} - v)- 

p 
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3.3.4 Complexity 0{n 2 ) of the algorithm 

As mentioned before, the complexity of this algorithm depends on the implementation of 
funcDescription data structure. 

The simplest implementation of funcDescription uses an array of pairs (v,p) sorted by 
v in ascending order and by p in the reverse order in the same time. This is possible based 
on Lemma 13. 3. 21 

The first block of instructions can be performed in G(# funcDescription) by simply 
shifting all elements such that p < to the right, and modifying v by iterating over all 
elements of funcDescription. The next block (computing of leftValue) can be computed 
with the same complexity, as min-{V : (v', _) S funcDescription A v' > v} is simply the 
next element after v in the ordered array. Finally, every instruction in the last block can 
be performed in ©(# funcDescription) by iterating over all its elements. 

In each iteration of the main loop at most 1 element is added to funcDescription. 
Therefore the computational complexity of the algorithm is 0(n 2 ) while the memory com- 
plexity is 0(n). 

3.3.5 Complexity O(nlogn) of the algorithm 

The previous result can be improved to 0(n log n) by using balanced binary search trees 
data structure. 

In this implementation funcDescription is represented by global variables Pmodifier and 
a BST of values (Av,p) where p is the key. An entry (v,p) in this structure is represented 
as 



E 



Av' ,p 

(Av' ,p')£funcDescriptionAp'^p 

so obtaining an element of funcDescription may take linear time. 

The advantages of this structure can be easily seen when analyzing the first block of 
the code. The division of funcDescription by the value of p (at first 0) can be achieved 
in O(logn). Shifting all elements of those subsets can then be done in a constant time 
by modifying first elements of those sets. Adding the extra node also requires O(logn) 
operations. 

Setting leftValue may require linear time, but all (apart from one) visited nodes will 
be removed in the third block. In all iterations of the main loop this instruction takes, 
therefore, 0(n). 

Removing nodes with the first coordinate ^ — 1 is obviously done in 0(n) in total. 
Identifying nodes with the first coordinate ^ 1 might seem problematic. It is, however, 
known that for the least p the value of v is equal to 1 + d. Relevant nodes can be, therefore, 
removed from the back in 0{n). Adding a n to the second coordinate of each node is done 
by adding it to global variable Pmodifier- 
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All iterations of the main loop require 0(n log n+n) operations. The memory complexity 
is also O(ralogn). 

3.3.6 Performance of the algorithm for the flat distance 

Performance of the algorithm depends on the choice of funcDescription data structure. 
Theoretic bounds for computational complexity are, however, not sufficient to argue about 
performance of these two options. The first reason is that the each operation in 0(n 2 ) 
algorithm is much faster than in 0(n log n) in terms of number of instructions. Secondly, 
hardware architectures provide solutions in which iterating over large tables is accelerated. 
Finally, the algorithm does reach its theoretical bound only if many points concentrate on 
a small interval. A gap of size 2 between two points completely cleans funcDescription 
data structure. All that is shown in numerical tests. To measure performance we have used 
a single core of AMD Athlon II X4 605e processor clocked at 2.3Ghz with 8GB of memory. 
The results are presented in Figsjl](2] 
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N Dirac deltas randomly distributed over [-1,1] 




20K 30K 



Figure 1: Comparison of the performance of the two proposed algorithms for the flat 
distance for N Dirac deltas randomly distributed over [—1,1]. The plot shows how the time 
of computation depends on N. For each input size 100 independent tests were executed 
to demonstrate how sensitive the algorithms are to input data distribution. Results of 
0(nlogn) algorithm are depicted as red dots, and results of 0(n 2 ) algorithm as blue dots. 
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N Dirac deltas (distances between each two: >2) 
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Figure 2: Comparison of the performance of the two proposed algorithms for the fiat 
distance for N Dirac deltas with a distributed over a large domain, i.e. distance between 
each two masses is larger than 2. In this case both algorithms are in fact linear, as the 
funcDescription structure has at most two elements. The plot demonstrates the overhead of 
using BST structures. Results of 0{nlogn) algorithm are depicted as red dots, and results 
of 0(n 2 ) algorithm as blue dots. 
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